Largest Lyapunov exponents for lattices of interacting classical spins 
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We investigate how generic the onset of chaos in interacting many-body classical systems is in the 
context of lattices of classical spins with nearest-neighbor anisotropic couplings. Seven large lattices 
in different spatial dimensions were considered. For each lattice, more than 2000 largest Lyapunov 
exponents for randomly sampled Hamiltonians were numerically computed. Our results strongly 
suggest the absence of integrable nearest-neighbor Hamiltonians for the infinite lattices except for 
the trivial Ising case. In the vicinity of the Ising case, the largest Lyapunov exponents exhibit 
a power-law growth, while further away they become rather weakly sensitive to the Hamiltonian 
anisotropy. We also provide an analytical derivation of these results. 
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The concept of microscopic chaos in many-particle sys- 
tems plays an essential role in the foundations of sta- 
tistical physics [ll-0|- In classical systems, chaos is de- 
fined by the appearance of exponential instabilities with 
respect to infinitesimal perturbations of the initial con- 
ditions. The spectrum of these instabilities is charac- 
terized by a set of eigenvalues - Lyapunov exponents - 
and the corresponding eigenvectors. In general, interact- 
ing, many-particle, classical systems are expected to be 
chaotic. The largest Lyapunov exponents and the whole 
Lyapunov spectra have been calculated numerically for 
classical many-body systems, such as gases of hard-core 
particles [1, , fluids with soft interactions Q , and lat- 
tice two-dimensional rotators 0, 0] , and analytically in a 
few cases |8hl4i]. Nevertheless, little is known about the 
universality of the features of the Lyapunov spectra, in 
particular when it comes to the systems with smooth dy- 
namics. It is also not clear how generic the onset of chaos 
actually is, and what happens in many-particle systems 
in the vicinity of integrable, i.e. nonchaotic, limits of the 
microscopic Hamiltonians. 

This Letter deals with the above issues on the basis 
of a systematic numerical and analytical investigation of 
the systems of interacting classical spins. These systems 
have been extensively studied, e.g., in the context of the 
spin diffusion problem isl - 17 1 , but their Lyapunov spec- 
tra have not yet been computed. Finite classical spin 
systems are also known to exhibit nontrivial integrable 
limits 18|. Another aspect of our motivation is the con- 
nection to the quantum case. There is a growing appre- 
ciation that generic interacting many-particle quantum 
systems exhibit relaxation behavior similar to their clas- 
sical chaotic counterparts, in particular as far as nuclear 
spin decays in solids are concerned 3-27|. In a broader 
context, the important conceptual development in the 
present study is the massive character of the numerical 
investigation covering an entire class of Hamiltonians. 

In this Letter, we consider several large lattices of in- 
teracting classical spins. For each lattice, we compute 



the largest Lyapunov exponent A,nax for several thousand 
randomly selected microscopic Hamiltonians. Thereby, 
we obtain the dependence of the largest Lyapunov expo- 
nent on the anisotropy of the spin-spin interaction. Our 
results give a strong indication of the absence of inte- 
grable Hamiltonians for infinite spin lattices with nearest- 
neighbor interaction besides the trivial Ising case (ex- 
plained below). We also find that, to the extent afforded 
by our numerical accuracy, the system becomes chaotic 
in the immediate vicinity of the integrable Ising case, 
with Amax exhibiting a universal power-law scaling. Fur- 
ther away from the Ising limit, A,nax becomes only weakly 
dependent on the Hamiltonian anisotropy, especially for 
bipartite lattices. The letter is concluded with a simple 
analytical derivation that describes the above-mentioned 
numerical results. 

We investigate the seven lattices shown in Fig.[TJ (LI) 
a chain, (L2) a rectangular ladder, (L3) a square lattice, 
(L4) a bilayer of square lattices, (L5) a cubic lattice, (L6) 
a triangular ladder and (L7) a triangular lattice. The 
interaction Hamiltonian for each lattice is of the nearest- 
neighbor (NN) type with periodic boundary conditions: 



NN 



H ^ ^ J X Six Sjx 



(1) 



where {Six, Siy, Siz) = Si are the three projections of the 
ith classical spin normalized by the condition Sf = 1. 

We numerically integrate the equations of motion as - 
sociated with the Hamiltonian ([T]): Si = Si x h; 
where hi is the local field given by the expression 



hi — ^ ^ JxSjx^x 4^ JyS'jy^y ^ J zS jz^z 
3{i) 



(2) 



Here e^;. By and are the unit vectors along the re- 
spective directions, and implies the summation over 
the nearest neighbors of the i-ih lattice site. We use the 
fourth-order Runge-Kutta algorithm and choose a time 
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FIG. 1. Lattices investigated in this work. Bipartite: (Ll) 
chain, (L2) rectangular ladder, (L3) square lattice, (L4) bi- 
layer of square lattices, (L5) cubic lattice. Non-bipartite: (L6) 
triangular ladder, (L7) triangular lattice. 



step of 0.005, sufficiently small so that, on the time scales 
of our simulations, energy is conserved. The initial con- 
ditions of each trajectory are chosen randomly on the 
energy shell with zero total energy, which corresponds to 
infinite temperature (29j . 

In order to obtain Amax, we numerically calculate a 
phase space trajectory and, at every time step, track the 
evolution of a tangent space vector that defines an in- 
finitesimal perturbation of it 22, [30| . The asymptotic 
average growth rate of this vector is equal to the largest 
Lyapunov exponent. In all cases, Amax becomes size- 
independent for sufficiently large lattices. 

For each lattice, we computed A,nax for many combi- 
nations of the coupling constants randomly selected on 
the "interaction sphere" J"^ -\- Jy -\- — 1. Specifically, 
8000 combinations selected from the isotropic distribu- 
tion were selected for the lattices (Ll), (L3), (L5), and 
2000 combinations for each of the remaining lattices. In 
addition, in order to investigate the scaling behavior near 
the Ising limit, we have processed more combinations in 
the vicinity of J2 = 1: 1000 for the lattices (Ll), (L3), 
and (L5) and 250 for the others. 

Our main numerical findings, namely A,nax as a func- 
tion of the parameter Jmax = Diax(| J^,!, | Jy|, | J^j), are 
presented in Fig. [2] The maximum value Jmax = 1 is 
realizable only in the Ising case, and thus represents the 
integrable limit with Amax = 0. The minimum value of 
Jmax is l/v^. It corresponds to either Heisenberg case 



Jx = Jy = Jz or "anti- Heisenberg" case Jx — Jy — — Jz 
(or equivalent cases). In Fig. [3l we present the rescaled 
plots Amax (Jmax)/ Amax (l/\/3) for lattices (L1-L5). On 
the basis of the results presented in Figs. [5] and |31 we 
make several important observations. 
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FIG. 2. (Color online) Largest Lyapunov exponents. Each 
point represents one Amax obtained numerically for a lattice 
indicated in the plot legend with one randomly chosen set 
of values J^, Jy and Jz as described in the text (a) Linear 
plot, (b) Log-log plot. The inset shows the prefactor a of 
the power-law fit Amax = a(l — Jmax)^''^ as a function of 
the number of nearest neighbors n with squares for (L1-L5), 
triangles for (L6,L7) and solid line for the fit a = n^''^. 



(i) For all lattices considered, no integrable cases be- 
sides the Ising case Jmax = 1 were found. This indi- 
cates that, for infinite lattices of classical spins with the 
nearest-neighbor interaction, the existence of another in- 
tegrable case is highly unlikely. 

(ii) The value of Amax is mainly controlled by Jmax, 
especially for the bipartite lattices. 

(iii) The dependence Amax(Jmax) for bipartite lattices 
(L1-L5) has nearly universal form, as illustrated by the 
rescaling shown in Fig. [S] 

(iv) The above dependence is nearly flat below Jmax ~ 
0.85; i.e., away from the integrable limit Amax is very 
weakly sensitive to the details of microscopic interaction. 

(v) Near the integrable limit Jmax = 1, to the best of 
our numerical accuracy (1 — Jmax > 10^^), each lattice, 
bipartite or not, becomes immediately chaotic, and, as 
shown in Fig. (H^b), exhibits an approximate power-law 
scaling Amax — a(l — Jmax)^^"^, whcrc a is a constant. 

(vi) As can be seen from Fig. HJa), the nonbipartite 
lattices (L6) and (L7) show a fork-shaped spread of Amax 
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FIG. 3. (Color online) Largest Lyapunov exponents for the 
bipartite lattices (L1-L5) from Fig. [2ja) rescaled by dividing 
by Ao = Amax(l/\/3). The dotted line represents Amax = 
Jmax(l — Jmax)^''*- The inset shows Ao as a function of the 
number of nearest neighbors n. 



as Jmax approaches 1 / \/3. The upper and the lower tips 
of the fork correspond to the anti-Heisenberg and Heisen- 
berg cases, respectively. 

Gross features of the above results can be reproduced 
by rather simple analytical estimates. In particular, away 
from the Ising case, the plateau values of Amax seen 
in Fig. ^a) at Jmax < 0.85 can be estimated with a 
factor-of-two accuracy as the typical frequency of one- 
spin motion given by the root-mean-squared value of 
the local field at the infinite temperature: Amax ^ 
+ Jy + j|)/3]^/^, where n is the number of the 
nearest neighbors. This expression also predicts that the 
heights of these plateaus scale as n^/^, while the inset 
of Fig. 13] indicates that the actual scaling is close but 
somewhat steeper. 

Now we turn to the approximation for the dependence 
Amax( Jmax) as the systcm approaches the Ising limit. We 
assume that Jz ^ Jx,Jy, which implies that Jmax = Jz ■ 
We also introduce variable Jj^ = [(Jj -f J^)/2]^/^ — 

[(1 — Jmax) /"A^^"^ denote the typical value of the trans- 
verse coupling. We consider two phase space trajectories 
{Si(<)} and {Si(i) + (5S,(t)}, where {5S,{i)} is an in- 
finitesimal difference. In order to estimate Amax, we lin- 
earize the equations of motion ([2]) with respect to small 
5Si{t) and keep the leading order in terms of J_lI Jz'- 



dt 

dSSi: 

dt 



= Jzj2cnt)ss,z , 

= Ji^C;^(<)<55,^ , 



(3) 
(4) 



where SStz and SSi^ are the projection of vectors SSi on 
the directions of vectors and x Si{t), respectively. 
The third projection of 6Si does not appear in Eqs. (jS]) 
and (HI), because it can be expressed in terms of SSiz - 



a consequence of the constraint Sf = 1. The parame- 
ters Cl^{t) and C^J{t) are determined by Si{t) and Sj{t) 
and have characteristic fluctuation times 1/J_l and l/Jz, 
respectively. 

Now, we make an assumption justified by the final re- 
sult [Eq. p^ ] that J± <^ Amax ^ Jz- We estimate the 
growth of the typical values of 6Siz and 6Siip over time 
T such that Amax ^ I/t ^ Jz- On the timescale r, 
the parameters Cl^ (t) stay nearly constant, while CJ^ {t) 
strongly fluctuate, so that {C1^{t))T ~ 0. We first write 



SS^z{t 



5S,^{t)+TJzY,Cl'{t)5S^ 

m 

ft + T 



(5) 



) = ss.ut) + J± fdt' C'^it')^S,^{t').{6) 



In this problem, a relatively slow growth of 5*^,^ is cou- 
pled to a random- walk- like growth of Siz- In order to 
extract Amax, we, therefore, look at the leading terms in 
the growth of SSf^ and SSf^: 

SSlit + r) « 6Sl{t) + 2TJzY.Cl^{t)SS,z{t)5S^^it),{7) 



SSlit + T) = SSfzit) + Jl 



3(i) 

t + T 



dt' 



t+T 



dt 



5Z ""v 



C:^{t')C^^{t')5S,^{t')5Sk^{t') - (8) 



In Eq.([8l), we neglected the term linear in J^, because 

(cj^-(t)).«o. 

Since we are only concerned with the scaling of Amax 
with Jmax and n, we convert Eqs. (O and (jH]) into 
an order-of-magnitude estimate for the typical growth 
of 5Sf^ and 5Sf^- The estimate includes (i) drop- 
ping lattice index in Eqs. ([7|) and ([8]), (ii) estimat- 
ing the instantaneous values of parameters C],^ (t) and 
C^{t) by 1, (ni) replacing the sum in Eq. ^ by ^/n 
(a consequence of the random sign of the n terms in 
that sum), and (iv) approximating the integral term in 
Eq. (ED by JlriSS^)^ dt'(E,(.) C;^'(OC^^ (i + t'))t ^ 
J\t{6S^Y \/nj J z- The factor ^JnjJz in the latter esti- 
mate is due to n terms in the sum, each producing a con- 
tribution to the integral of the order Xjhiz ^ l/{Jz-Jn). 
Finally, after dividing thus simplified Eqs. ([7]) and ([8]) by 
{SStp)^ and [SSz)^, respectively, we obtain 



5Sl{t) 

5Sl{t + T) 
5 Slit) 



1 -I- 2T^/nJ, 



1 



SSzjt) 

5s^{ty 

-Jj 5 Slit) 
\jz 5 Slit) 



(9) 



(10) 



Since the parameters dS^, and 5 Si should grow at the 
same rate, this implies that 

5Szit) ^ ^ ^ 



SS^it) 



(1 </max) 



1/3 



(11) 
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FIG. 4. Projections of the Lyapunov vector correponding to 
Amax onto the subspaces {SSiz} and {SSi^p} for the 16 x 16 x 16 
cubic lattice (L5). Solid line is the fit based on Ea.([TT|l. 



Substituting Eq. ([TT]) into Eq. ^ and comparing the 
right-hand side with the expression l + 2Amax''', we finally 
obtain the estimate: 



Amax — "^max "A 



Nl/3 
-'max/ 



/^(l-Jmax)'/', (12) 



which is consistent with the numerically observed 1/3- 
power law shown in Fig.[5Jb). In Ref. [l3| the same power 
law was obtained for weakly interacting dilute gases via a 
perturbation expansion around the integrable ideal gas. 
The prefactor scaling as -^/n is also compatible with the 
numerical results, as illustrated in the inset of Fig.[2fb). 

Equation pT|) predicts further that the components 
of the Lyapunov eigenvector corresponding to Amax de- 
pend systematically on whether they are parallel to the 
z-direction or not. As illustrated in Fig. 01 this prediction 
also agrees with our numerical results. 

Our derivation of Eq. (fT2)) is based on the assumption 
of the random fluctuations of the sum in Eq. ^ . How- 
ever, in the Ising limit, the Fourier transform of this sum 
contains only a finite number of frequencies. Therefore, 
in the vicinity of the Ising limit, recurrences may occur 
if correlations in Siz do not decay sufficiently fast. Such 
recurrences, in turn, would contradict our assumption of 
the fast decay of the time correlations of the above sum. 
This would be most problematic for lattices with smaller 
numbers of nearest neighbors, such as spin chains. In- 
deed, we observe in Fig. [2][b), that spin chains exhibit 
larger deviations from the 1 /3-power law than other lat- 
tices. The above recurrences may, in fact, be the route to 
breaking down the chaotic nature of the spin dynamics 
around the Ising limit. However, our numerical results 
indicate that, if such a breakdown occurs, it occurs at 
values of \J±/Jz\ < 10"''. 

Finally, we mention that the above estimate may be 
repeated for the case of Jx and Jy smaller but not much 
smaller than J^. In this a case, r must be chosen much 
shorter than 1 /J± and thus both Cl^ (t) and CJ* (t) can be 



assumed constant on the time scale of r. This estimate 
would then give Amax = (nJmax)^/^(l - -'^^ax)^'^'': which, 
as illustrated in Fig. [31 exhibits a good overall agreement 
with the numerical results for bipartite lattices over the 
entire range of Jmax- 

In conclusion, we have presented a systematic study of 
largest Lyapunov exponents Amax for a very large variety 
of classical spin systems. Our findings strongly suggest 
the absence of integrable nearest-neighbor Hamiltonians 
for the type of lattices considered except for the Ising 
case. As far as the behavior of Amax is concerned, a 
number of universal features enumerated above as (ii)- 
(vi) are observed. We have also analytically derived the 
scaling of Amax with the anisotropy of the Hamiltonian. 
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